home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Cream of the Crop 20
/
Cream of the Crop 20 (Terry Blount) (1996).iso
/
os2
/
sysb091a.zip
/
sysbench
/
src
/
pmb_linpack.c
< prev
next >
Wrap
C/C++ Source or Header
|
1994-11-05
|
27KB
|
1,232 lines
/*
Translated to C by Bonnie Toy 5/88
You MUST specify one of -DSP or -DDP to compile correctly.
You MUST specify one of -DROLL or -DUNROLL to compile correctly.
You MUST specify a timer option(see below) to compile correctly.
To compile double precision version for Sun-4:
cc -DUNIX -DDP -DROLL -O4 clinpack.c
To compile single precision version for Sun-4:
cc -DUNIX -DSP -DROLL -O4 -fsingle -fsingle2 clinpack.c
To obtain rolled source BLAS, add -DROLL to the command lines.
To obtain unrolled source BLAS, add -DUNROLL to the command lines.
PLEASE NOTE: You can also just 'uncomment' one of the options below.
*/
/* #define SP */
#define DP
#define ROLL
/* #define UNROLL */
/***************************************************************/
/* Timer options. You MUST uncomment one of the options below */
/* or compile, for example, with the '-DUNIX' option. */
/***************************************************************/
/* #define Amiga */
/* #define UNIX */
/* #define UNIX_Old */
/* #define VMS */
/* #define BORLAND_C */
/* #define MSC */
/* #define MAC */
/* #define IPSC */
/* #define FORTRAN_SEC */
/* #define GTODay */
/* #define CTimer */
/* #define UXPM */
#ifdef SP
#define REAL float
#define ZERO 0.0
#define ONE 1.0
#define PREC "Single "
#endif
#ifdef DP
#define REAL double
#define ZERO 0.0e0
#define ONE 1.0e0
#define PREC "Double "
#endif
#define NTIMES 100
#ifdef ROLL
#define ROLLING "Rolled "
#endif
#ifdef UNROLL
#define ROLLING "Unrolled "
#endif
//#include <stdio.h>
#include <math.h>
static double st[8][6];
double pmb_linpack ()
{
static REAL aa[200][200],a[200][201],b[200],x[200];
REAL cray,ops,total,norma,normx;
REAL resid,residn,eps;
REAL epslon(),kf;
double t1,tm,tm2,dtime();
static int ipvt[200],n,i,ntimes,info,lda,ldaa,kflops;
lda = 201;
ldaa = 200;
cray = .056;
n = 100;
/*
fprintf(stdout,ROLLING);fprintf(stdout,PREC);
fprintf(stdout,"Precision Linpack\n\n");
*/
// fprintf(stderr,ROLLING);fprintf(stderr,PREC);
// fprintf(stderr,"Precision Linpack\n\n");
ops = (2.0e0*(n*n*n))/3.0 + 2.0*(n*n);
matgen(a,lda,n,b,&norma);
t1 = dtime();
dgefa(a,lda,n,ipvt,&info);
st[0][0] = dtime() - t1;
t1 = dtime();
dgesl(a,lda,n,ipvt,b,0);
st[1][0] = dtime() - t1;
total = st[0][0] + st[1][0];
/* compute a residual to verify results. */
for (i = 0; i < n; i++) {
x[i] = b[i];
}
matgen(a,lda,n,b,&norma);
for (i = 0; i < n; i++) {
b[i] = -b[i];
}
dmxpy(n,b,n,lda,x,a);
resid = 0.0;
normx = 0.0;
for (i = 0; i < n; i++) {
resid = (resid > fabs((double)b[i]))
? resid : fabs((double)b[i]);
normx = (normx > fabs((double)x[i]))
? normx : fabs((double)x[i]);
}
eps = epslon((REAL)ONE);
residn = resid/( n*norma*normx*eps );
/* printf(" norm. resid resid machep");
printf(" x[0]-1 x[n-1]-1\n");
printf("%8.1f %16.8e%16.8e%16.8e%16.8e\n",
(double)residn, (double)resid, (double)eps,
(double)x[0]-1, (double)x[n-1]-1);
fprintf(stderr," times are reported for matrices of order %5d\n",n);
fprintf(stderr," dgefa dgesl total kflops unit");
fprintf(stderr," ratio\n");
*/
st[2][0] = total;
st[3][0] = ops/(1.0e3*total);
st[4][0] = 2.0e3/st[3][0];
st[5][0] = total/cray;
// fprintf(stderr," times for array with leading dimension of%5d\n",lda);
print_time(0);
matgen(a,lda,n,b,&norma);
t1 = dtime();
dgefa(a,lda,n,ipvt,&info);
st[0][1] = dtime() - t1;
t1 = dtime();
dgesl(a,lda,n,ipvt,b,0);
st[1][1] = dtime() - t1;
total = st[0][1] + st[1][1];
st[2][1] = total;
st[3][1] = ops/(1.0e3*total);
st[4][1] = 2.0e3/st[3][1];
st[5][1] = total/cray;
matgen(a,lda,n,b,&norma);
t1 = dtime();
dgefa(a,lda,n,ipvt,&info);
st[0][2] = dtime() - t1;
t1 = dtime();
dgesl(a,lda,n,ipvt,b,0);
st[1][2] = dtime() - t1;
total = st[0][2] + st[1][2];
st[2][2] = total;
st[3][2] = ops/(1.0e3*total);
st[4][2] = 2.0e3/st[3][2];
st[5][2] = total/cray;
ntimes = NTIMES;
tm2 = 0.0;
t1 = dtime();
for (i = 0; i < ntimes; i++) {
tm = dtime();
matgen(a,lda,n,b,&norma);
tm2 = tm2 + dtime() - tm;
dgefa(a,lda,n,ipvt,&info);
}
st[0][3] = (dtime() - t1 - tm2)/ntimes;
t1 = dtime();
for (i = 0; i < ntimes; i++) {
dgesl(a,lda,n,ipvt,b,0);
}
st[1][3] = (dtime() - t1)/ntimes;
total = st[0][3] + st[1][3];
st[2][3] = total;
st[3][3] = ops/(1.0e3*total);
st[4][3] = 2.0e3/st[3][3];
st[5][3] = total/cray;
print_time(1);
print_time(2);
print_time(3);
matgen(aa,ldaa,n,b,&norma);
t1 = dtime();
dgefa(aa,ldaa,n,ipvt,&info);
st[0][4] = dtime() - t1;
t1 = dtime();
dgesl(aa,ldaa,n,ipvt,b,0);
st[1][4] = dtime() - t1;
total = st[0][4] + st[1][4];
st[2][4] = total;
st[3][4] = ops/(1.0e3*total);
st[4][4] = 2.0e3/st[3][4];
st[5][4] = total/cray;
matgen(aa,ldaa,n,b,&norma);
t1 = dtime();
dgefa(aa,ldaa,n,ipvt,&info);
st[0][5] = dtime() - t1;
t1 = dtime();
dgesl(aa,ldaa,n,ipvt,b,0);
st[1][5] = dtime() - t1;
total = st[0][5] + st[1][5];
st[2][5] = total;
st[3][5] = ops/(1.0e3*total);
st[4][5] = 2.0e3/st[3][5];
st[5][5] = total/cray;
matgen(aa,ldaa,n,b,&norma);
t1 = dtime();
dgefa(aa,ldaa,n,ipvt,&info);
st[0][6] = dtime() - t1;
t1 = dtime();
dgesl(aa,ldaa,n,ipvt,b,0);
st[1][6] = dtime() - t1;
total = st[0][6] + st[1][6];
st[2][6] = total;
st[3][6] = ops/(1.0e3*total);
st[4][6] = 2.0e3/st[3][6];
st[5][6] = total/cray;
ntimes = NTIMES;
tm2 = 0;
t1 = dtime();
for (i = 0; i < ntimes; i++) {
tm = dtime();
matgen(aa,ldaa,n,b,&norma);
tm2 = tm2 + dtime() - tm;
dgefa(aa,ldaa,n,ipvt,&info);
}
st[0][7] = (dtime() - t1 - tm2)/ntimes;
t1 = dtime();
for (i = 0; i < ntimes; i++) {
dgesl(aa,ldaa,n,ipvt,b,0);
}
st[1][7] = (dtime() - t1)/ntimes;
total = st[0][7] + st[1][7];
st[2][7] = total;
st[3][7] = ops/(1.0e3*total);
st[4][7] = 2.0e3/st[3][7];
st[5][7] = total/cray;
/* the following code sequence implements the semantics of
the Fortran intrinsics "nint(min(st[3][3],st[3][7]))" */
/*
kf = (st[3][3] < st[3][7]) ? st[3][3] : st[3][7];
kf = (kf > ZERO) ? (kf + .5) : (kf - .5);
if (fabs((double)kf) < ONE)
kflops = 0;
else {
kflops = floor(fabs((double)kf));
if (kf < ZERO) kflops = -kflops;
}
*/
if ( st[3][3] < ZERO ) st[3][3] = ZERO;
if ( st[3][7] < ZERO ) st[3][7] = ZERO;
kf = st[3][3];
if ( st[3][7] < st[3][3] ) kf = st[3][7];
kflops = (int)(kf + 0.5);
// fprintf(stderr," times for array with leading dimension of%4d\n",ldaa);
print_time(4);
print_time(5);
print_time(6);
print_time(7);
// fprintf(stderr,ROLLING);fprintf(stderr,PREC);
// fprintf(stderr," Precision %5d Kflops ; %d Reps \n",kflops,NTIMES);
return kflops;
}
/*----------------------*/
static print_time (row)
int row;
{
/*fprintf(stderr,"%11.2f%11.2f%11.2f%11.0f%11.2f%11.2f\n",
(double)st[0][row], (double)st[1][row], (double)st[2][row],
(double)st[3][row], (double)st[4][row], (double)st[5][row]);
*/
}
/*----------------------*/
static matgen(a,lda,n,b,norma)
REAL a[],b[],*norma;
int lda, n;
/* We would like to declare a[][lda], but c does not allow it. In this
function, references to a[i][j] are written a[lda*i+j]. */
{
int init, i, j;
init = 1325;
*norma = 0.0;
for (j = 0; j < n; j++) {
for (i = 0; i < n; i++) {
init = 3125*init % 65536;
a[lda*j+i] = (init - 32768.0)/16384.0;
*norma = (a[lda*j+i] > *norma) ? a[lda*j+i] : *norma;
}
}
for (i = 0; i < n; i++) {
b[i] = 0.0;
}
for (j = 0; j < n; j++) {
for (i = 0; i < n; i++) {
b[i] = b[i] + a[lda*j+i];
}
}
}
/*----------------------*/
static dgefa(a,lda,n,ipvt,info)
REAL a[];
int lda,n,ipvt[],*info;
/* We would like to declare a[][lda], but c does not allow it. In this
function, references to a[i][j] are written a[lda*i+j].
*/
/*
dgefa factors a double precision matrix by gaussia